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Abstract 

We propose a variety of models of random walk, discrete in space and time, 
suitable for simulating stable random variables of arbitrary index a (0 < a < 2), in 
the symmetric case. We show that by properly scaled transition to vanishing space 
and time steps our random walk models converge to the corresponding continuous 
Markovian stochastic processes, that we refer to as Levy-Feller diffusion processes. 
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1 Introduction 



By a Levy-Feller diffusion process we mean a Markovian process governed by 
a stable probability density function (pdf) evolving in time, g a {x, t; 9) , whose 
spatial Fourier transform (the characteristic function) reads 

+00 

0) = Je tKX g a (x,t;9)dx = exp (-t|*f e^ 11 «0^/2) , (1.1) 

—00 

where x , k £ R , t > . The two relevant parameters, a , called the index 
of stability, and 9 (related to the asymmetry) , improperly referred to as the 
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skewness, are real numbers subject to the conditions, see e.g. [|TJ, 

n ^ o \n\ ^ { a i if < a < 1 , n _v 

0<a ^ 2 ' '^U-a, if l<a<2. (L2) 

By introducing the similarity variable we can write g a (x,t;9) = 

t* 1 /" pafet* 1 / ', 9) , where p a (x;9) is the stable pcff at i = 1 . The specific 
form of the characteristic function (1.1) allows us to recognize g a (x,t;9) as 
the Green function (fundamental solution) of the Cauchy problem 

d 

— u(x,t) = Dg[u(x,t)], u(x,0)=5(x), xGR, £>0, (1.3) 

where Dg is the pseudo-differential operator with symbol 

^« = _| K |a e i(signK)07r/2_ (L4) 

Let us recall that a generic pseudo-differential operator A, acting with re- 
spect to the variable x G R, is defined through its Fourier representation, 
namely f^™ e lKX A [(j>(x)] dx = A(k) <P{k) , where <j>{x) denotes a sufficiently 
well-behaved function in R, and A(k) is referred to as symbol of A , given as 
A{k) = (Ae~ iKX ) e +iKX . 

With the names of Levy and Feller we have intended to honour both Paul 
Levy, M, who first introduced the class of stable distributions, see M, M, 
and William Feller 0, who first investigated the semigroups generated by a 
pseudo-differential equation of type (1.3-4). For a = 2 and a = 1 (with 9 = 0) 
we recover the standard Gaussian and Cauchy pdf's 

^(x,t;0) = -^t- 1/2 exp(-^-) , g 1 ( Xl t;0) = -^—^. (1.5) 

2^71 \ At J 7T X Z + t A 

In physics, the first recognition that the (symmetric) Levy distribution could 
be characterized via a pseudo-differential operator of type (1.3-4) was made ex- 
plicitly by West & Seshadri 0. Recently, Gorenflo & Mainardi 0, [§], 0, have 
revised Feller's original arguments by interpreting (1.3) as a space-fractional 
diffusion equation (of order a and " skewness" 9) and have provided a variety 
of related random walk models, discrete in space and time, which by properly 
scaled transition to vanishing space and time steps converge to the correspond- 
ing continuous Levy-Feller processes. In other words the discrete probability 
distributions generated by the random walk models have been proved to be- 
long to the domain of attraction of the corresponding stable distribution. 

Here, limiting ourselves to the symmetric case (9 = 0), we present the main 
features of the random walk models by Gorenflo & Mainardi, and we display 
preliminary results of a few numerical case studies, which can be of some inter- 
est in econophysics. The Levy statistics in modelling fluctuations of economical 
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and financial variables, formerly used by Mandelbrot in the early sixties, is 
still nowadays adopted with success, see e.g. flO ], [fTTfl , [ 12| . 



2 Outline of the general theory 

Let Y be an integer-valued random variable and let the random variables 
Yi,Y 2 ,Y 3 , ... be i.i.d. (= independent identically distributed), all having their 
probability distribution common with Y . We define a spatial-temporal grid 
{(xj , t n ) \j G Z, n G N } by = = jh, t n = t n {r) = nr , where 

h > and r > . Then we consider the sequence of random variables 

S n = hY 1 + hY 2 + ... + hY n , n G N , (2.1) 

with (for convenience) Sq = , and interpret it as follows. A particle, sitting 
in x = xq = at time t = to = finds itself at a later instant t — t n in point 
x — S n which is an integer multiple of ft . We recognize the pu = P(Y = k) 
(for k G Z) as "transition probabilities": p& is the probability of a particle 
jumping from a point Xj — S n to a point = 5* n+ i as time proceeds from 
t n to t n+ i . All p fc are non-negative, and their sum equals 1. 

The probability yj(t n ) of sojourn of our particle in point Xj at instant t n obeys 
the transition law 

+oo 

%(Wi) = X) PkV3-k{tn): %'(0) = 5 J o, jGZ, Ii6N . (2.2) 

fc=— oo 

which has the form of a discrete convolution. Hence, introducing generating 
functions 

P( Z ) = Pi* 1 ' Vn( Z ) = Vjitn)^, (2-3) 
j=— oo j=—oo 

we obtain 

&(z) = yo(z) ■ [p(z)} n = [p(z)] n , n G N . (2.4) 

The power series in (2.3) and (2.4) are absolutely and uniformly convergent on 
\z\ — 1 and assume the value 1 at z = 1 . Putting z = e lK , /t G R , and observ- 
ing = e lK i h = e lKXj , we recognize p(n; h) = p(e lKh ) and y(n, t n ; h) = y n (e %Kh ) 
as characteristic functions of the random variables h Y and S n , respectively. 

Our aim is to approximate the Levy-Feller diffusion process, which is governed 
by the evolution equation (1.3), arbitrarily well. To this purpose we introduce 
a strictly monotonic scaling relation r = o~{h) — ► as h — > . We will fix 
t > and let h (and likewise r) go to zero over such values that always 
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n = t/r = t/a(h) is a positive integer. Then we have the equivalences 



n — > oo •<=>- ft — > •<==>- r — > , 

and ft, depends on r, finally on n, so that ft = h(n) . Replacing ft by h{n) in (2.1) 
we obtain a sequence of random variables X n with characteristic functions 
y{K,t n ;h) = [p{e iKh )} n (note that now t n = t is fixed). Invoking Theorem 3.6.1 
of Lukacs [O], what remains to be shown is that y{n,t; ft) — > exp(— t|fi;| a ) as 
ft — > , the characteristic function of the corresponding symmetric Levy-Feller 
process. For this it suffices that, for fixed k^O, 

\og[y( K ,t;h)] = -^-log\p(e iKh )]^-t\K\ a , as ft^O. (2.5) 
a(h) 

Our random walk can be interpreted as a "difference scheme" to approximate 
the evolution equation( 1.3), if we write (2.2) in the equivalent form, observing 
the scaling relation, 

yj (tn + l) ~ Vjitn) _ 1 [( po -l)y j (t n )+^ Pk y j _ k (t n )] 1 (2.6) 



a(h) 



In fact, whereas the L.H.S. is the explicit discrete approximation to the first- 
order time derivative jLu{x, t) , the R.H.S can be considered as a particular 
discrete approximation to the space pseudo-differential term Dq [u(x, t)] , pro- 



vided we mean Uj(t n ) = Jx--h/2 u ( x ^n) dx ~ hu(xj,t n ) with Uj(0) = 5 



jo 



3 The random walk models 



From the previous Section we have learnt that, in order to construct discrete 
random walk models which are convergent (in distribution) to the (symmetric) 
stable pdf's, the clue points are: 1) to guess a suitable generating function 
p(z) , whose coefficients of its power series expansion provide the transition 
probabilities, 2) to determine the corresponding scaling relation r = cr(ft) 
which ensures the required convergence. 

In the classical case of the Gaussian distribution (a = 2) the matter is eas- 
ily treated if we remember that the corresponding density is the fundamental 
solution of the standard diffusion equation, which is known to be well approx- 
imated via the finite-difference equation 

Vj{t n+1 ) - Vjitn) y j+ iit n )-2y j {t n )+y j _ 1 (t n ) 

7a ' %l u J- d io- (3.1J 

In this case, introducing the scaling parameter \i = r/h 2 , so r = cr(ft) = \i ft 2 , 
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the transition probabilities turn out to be 



p = 1 - 2/x , p±i = (i , p± k = , k = 2, 3 . . . , . 



(3.2) 



subject to the condition < /x < 1/2 . Thus the generating function is 



The proof of the convergence to the Gaussian is simple since one easily finds 
[t/(nh 2 )]log[p(e tKh ] — > —tn 2 as h — » . The scheme (3.2) means that for 
approximation of the standard Gaussian process the corresponding random 
walk model exhibits only jumps of one step to the right or one to the left or 
jumps of width zero. For the stable non Gaussian processes we expect to find a 
non-polynomial generating function with infinitely many transition coefficients 
which imply the occurrence of arbitrarily large jumps. It is common practice 
to refer to the corresponding random walks as to Levy flights. 

In the following, limiting ourselves to the symmetric cases (8 = 0), we shall 
resume the main features of three different random walk models, referred to 
as (RW1), (RW2), (RW3), of which Gorenflo and Mainardi have proved the 
convergence to the corresponding continuous processes. For each model we 
give the generating function p(z) with the transition probabilities pk and the 
scaling relation r = a(h), referring to the original papers for details. From the 
analysis of the classical Gaussian case we find it natural to first introduce the 
scaling parameter /x = r/h a with < a < 2 but, as we shall show later, this 
will not necessarily imply u(h) = /x h a , with /x constant for fixed a . 

The model (RW1) has been introduced and discussed in J7j, starting from 
the identification of the operator Dq in the framework of fractional calculus, 
see ||, and then applying, in the authors' original approach, the Griinwald- 
Letnikov discretized scheme. For this model we need to keep distinct the two 
cases (a) < a < 1 and (b) 1 < a < 2 , being the case a = 1 excluded in this 
treatment. The generating function is 



For both cases the scaling relation is confirmed to be r = a{h) = /x h a , but 
the parameter /x is subject to different restrictions to ensure that < p < 1 . 

In the case (a) we have < /x < cos(o;7r/2) and 



p(z) = l + fi[z-2 + z 1 }. 



(3.3) 




< a < 1 , 



1 < a < 2. 



(3.4) 




for k 



1,2,3,... . 



(3.5a) 
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In the case (b) we have < /i < | cos(a7r/2)|/a and 



PO = 1 - 1 7 TTyTT > P±l ~ 




cos(a7r/2)| ' 2| cos(a7r/2)| 

p±k = (-l)k+i ^ - - ( i a ) , for A; = 2,3,4,... 

^ V 7 2|cos(a7r/2)| U+ V 



(3.56) 



We note that, whereas the classical Gaussian random walk (3.2) is promptly 
recovered from (3.5b) for a = 2, the random walk for the Cauchy process 
(a = 1) cannot be obtained, neither directly nor by a passage to the limit 
a — *■ 1 . Indeed, in both the limits a — > 1 _ and a — ► 1 + the permissible range 
of the scaling factor \i is vanishing. In numerical practice the consequence will 
be that if a is near 1 the convergence is slow: for good approximation we will 
need a very small step-time r with respect to the step-length h . 

Differently from (RW1) we shall show how the two other random walk models, 
(RW2) and (RW3), that we are going to briefly discuss, exhibit a smooth 
scaling law in the range < a < 2, so the case a = 1 is no longer singular. 
However, "there is no free lunch". We have to pay for the good behaviour at 
a = 1 with bad behaviour at a = 2 in a sense to be seen later. 



The model (RW2) is obtained by imposing for the transition probabilities an 
expression in terms of binomial coefficients, as suggested from (3.5b), which 
is required to be regular in the limit as a — > 1 , namely 

Po = l-2A, p±* = (jfe+l)' i ovk = 1,2,3,... , (3.6) 

where A , subject to the condition < A < 1/2 is to be determined. The model 
has been extensively discussed in whereas the particular case a = 1 (related 
to the Cauchy process) has been formerly presented in ||. To guarantee the 
convergence for all a (0 < a < 2) to the corresponding continuous process, 
the asymptotics as h — > still requires a scaling relation of the kind r = \ih a 
(with \i constant), namely 

f 2Acos(mr/2) , 
i — — h a if 0<a<2, a^l, 

I A h 2 if a = 2 . 

We recognize that for < a < 2 the model (RW2) exhibits a smooth scaling 
law but a discontinuity is present at a = 2 as can be seen by taking the limit 
as a — >• 2~ . This allows us to recover for a = 2 the Gaussian model (3.2)-(3.3), 
but it shows that in numerical practice when a is near 2 the convergence is 
expected to be slow. For < a < 2 we obtain the generating functions and 
the transition probabilities as follow. Putting p(z) = [(1 — z) a ~ l — 1] for a ^ 1 , 
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and introducing the scaling parameter \i , we have 

p(«) = 2 cos(a7r/2) (3 _ g) 

1 - - [(1 - z- 1 ) log (1 - z) + (1 - *) log (1 - z-% a = l. 

(. IT 

For < a < 2 , 1 we have < /x < cos(a7r/2)/(l — a) and 

— a) 



Po = 1 



cos(«7r/2) 



p ±fc = (-l) fc if [ a I for fc = 1,2,3,... 

F±fc V ; 2cos(avr/2) \fc + 1/ ' ' ' ' 



(3.9) 



for a = 1 we have < fi < it/ 2 and 
2/i fi 1 

7T ' ^ ±fc 7T k(k + 1) 



Po = l~ — , P±k = - tjt~. 7s , for fc = 1,2,3,4,... . (3.10) 



In both models (RW1), (RW2) the generating function is expressed in terms 
of elementary functions and the transition coefficients, for < a < 2 , exhibit 
an asymptotic behaviour as \k\ — > oo consistent with that of the power- law 
tails of the stable densities as \x\ — > oo , see e.g. [|TJ. Indeed we obtain 

p k ~/ir(a+ 1) sm ( 7ra / 2 ) |fc|-(°+i) as — , oo , 0<a<2; (3.11) 

7T 

Pa {x; 0) ~ T(a + 1) g^i!^ZS |a;|-( a+1 > as |x| -»• oo , 0<a<2. (3.12) 

7T 

The asymptotic behaviour of the stable densities is the starting point for the 
random walk model (RW3) in that, following Gillis & Weiss [jl~4| , we require 
that all pk for k ^ are proportional to |A;|~( Q!+1 ' ) . However, the simplicity 
of the starting point leads to difficulties for treating this model since the 
generating function is no longer elementary, a fact that makes the convergence 
proof a really hard affair. We need to recall the following special functions 

OO OO 

C(0) = (3>1, 9(z,0) = y £^g, \z\<l, PER, (3.13) 

fc=i k=i K 

respectively known as the Riemann zeta function and the polylogarithmic func- 
tion. This model has been extensively discussed in M (see also WM), where 
the following generating function is derived 

p(z) = 1 - 2A C(a + 1) + A a + 1) + a + 1)] . (3.14) 

Here = a + 1 > 1 so $(2,/?) is by its power series also defined on the 
periphery \z\ = 1 of the unit circle. Then we get a pure power-law random 
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walk with 



Po = l-2AC(a+l), Pk = Mk\ 



-(a+l) 



for k ^ . 



(3.15) 



where A, subject to the condition < A < l/[2£(a + 1)] , is to be deter- 
mined. To ensure the convergence to the corresponding continuous process 
the following scaling relation must hold 



a{h) 



Att 



T(a + 1) sin(a7r/2) 
\h 2 \\ogh\ 



if 
if 



< a < 2 



(3.16) 



a 



2. 



It is interesting to note that in the case a = 2 the classical random walk model 
(3.2)-(3.3) is no longer recovered, since now arbitrarily large jumps occur (with 
a probability decaying as k~ 3 ) as in the Levy flights. Nevertheless, through the 
scaling relation (3.16), this random walk converges to the continuous Gaussian 
process. Thus this model gives us the opportunity to verify that a random 
walk with infinite variance steps may (slowly) converge to the Gaussian, in 
agreement with a general theorem on the domain of attraction of the normal 
law, see e.g. fl5|. 



4 Numerical results 



In general the random walk models are not only valuable from the conceptual 
point of view for visualizing what the diffusion means but also for numerical 
calculations, either as Monte Carlo simulation of particle paths in a diffusion 
process or as discrete imitation of the process in form of redistribution (from 
one time level to the next) of clumps of an extensive quantity (across the 
spatial grid points). Our models can be used in at least three different ways: 
(a) as finite difference schemes for approximate calculation of symmetric stable 
densities; (b) for producing sample paths of individual particles performing the 
random walk; (c) for producing histograms of the approximate realization of 
the densities g a by simulating many individual paths with the same number 
of time steps and making statistics of the final positions of the particles. 



For numerical simulations of stable random variables different algorithms have 
been provided by a number of specialists, including Chambers et al ||16|| , Bar- 
tels Mantegna fL8| , Janicki & Weron |19| , Samorodnitsky & Taqqu |20 



Our present approach for treating Levy statistics has been carried out inde- 
pendently from all the above references but uniquely based on the random 
walk models presented here, so, as far as we know, our results would be in the 
great part original. 

Having preliminarily checked a sufficient level of accuracy for our finite differ- 
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200 400 600 800 1000 1200 -3 -2 -1 1 2 3 

time steps x 

Fig. 1. Sample path (left) and histogram (right) for a = 1 . (Cauchy) 



a -1.5, final time - 1.33 a -1.5 




time steps x 



Fig. 2. Sample path (left) and histogram (right) for a = 1.5. 

ence schemes with the existing tables of stable densities, here we present some 
results on the simulation of the sample paths and histograms corresponding 
to some typical values of the index of stability, namely a — 1, 1.5, 2 . In prac- 
tice, in our numerical studies there is required truncation in two ways. It is 
impossible to simulate all infinitely many discrete probabilities, so the size 
of possible jumps must be limited to a maximal possible jump length. The 
other truncation is required if a priori one wants a definite region of space 
to be considered in which the walk takes place. Then, particles leaving this 
space have been ignored. Our simulations, based on one million of realiza- 
tions, have been carried out in the interval \x\ < 4 . All the histograms refer 
to stable densities at t — 1 for \x\ < 3 , the space interval being reduced to 
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time steps 



Fig. 3. Sample path (left) and histogram (right) for a = 2 . (Gauss) 

avoid the border effects. The sample paths are plotted against the time steps, 
up to 1200 for a = 1 and up to 2000 for a = 1.5 , 2 , so they refer to different 
final times, namely t = 1.87, 1.33,0.8, respectively. The transition probabili- 
ties have been chosen from our random walk models as follows: a = 1 from 
(RW2), see (3.10), with scaling parameter fi = 7r/4; a = 1.5 , from (RW1), 
see (3.5b), with \i = (2/3) cos(37r/4) ; for a = 2 , from the standard model, see 
(3.2), with /i = 1/4. The cases a = 1 (Cauchy process) and a = 2 (normal 
process) have been considered for a possible comparison with the standard 
and accurate algorithms existing in the literature, whereas a = 1.5 has been 
chosen in view of possible applications in econophysics where usually the index 
of stability ranges from 1.4 to 1.7, see e.g. pT0|] , [ PH . 



5 Conclusions 

For the simulation of Markovian processes characterized by symmetric Levy 
probability densities evolving in time we have presented three different random 
walk models, discrete in space and time, by giving their respective transition 
probabilities. We have indicated how, by use of generating functions, conver- 
gence for properly scaled transition to vanishing steps of space and time can be 
analyzed, and we have hinted at peculiarities occurring at the particular sta- 
bility indices a = 1 (the Cauchy process) and a = 2 (the Gauss process). We 
have displayed preliminary results of a few numerical case studies concerning 
sample paths and histograms to check the efficiency of our algorithms. From 
the sample paths one can recognize the "wild" character of the Levy flights 
with respect to the "tame" character of the Brownian motion. 
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We expect that our arguments can be relevant in different fields of physics 
including the emerging one of econophysics, where stable distributions are be- 
coming more common. In statistical physics the stable distributions play a key 
role in the (wonderful) world of random walks constructed by the late Mon- 
troll and continued through his school, see e.g. [^TJ, [p^j , [ f23fl . Also Tsallis and 
his associates have recognized the key role of stable distributions with respect 
to a generalized theory of thermostatics, see e.g. |24||, and references therein. 



Here, we have (only briefly) pointed out the relation between Levy statistics 
and space-fractional diffusion equations. However, stable distributions turn 



out to be related also with time-fractional diffusion equations, see e.g. p5 
Furthermore, this topic is relevant for fractal phenomena, where differential 
equations of fractional order are usually adopted to describe their evolution, 



see e.g. [26|, | 27fl , [28fl, and reference herein. 



Acknowledgements 



We gratefully acknowledge fruitful discussions with Rosario Mantegna. This 
work was supported in part by the Italian CNR and INFN, and by the Research 
Commission of Free University of Berlin. 



References 



[1] W. Feller, An Introduction to Probability Theory and its Applications, Vol. 2 
(Wiley, New York, 1971). 

[2] P. Levy, Theorie des erreurs. La lois de Gauss et les lois exceptionnelles, Bull. 
Soc. Math. France 52 (1924) 49-84. 

[3] P. Levy, Calcul des probabilities (Gauthier-Villars, Paris, 1925). 

[4] P. Levy, Theorie de V Addition des Variables Aleatoires (Gauthier-Villars, Paris, 
1937). 

[5] W. Feller, On a generalization of Marcel Riesz' potentials and the semi-groups 
generated by them, Meddelanden Lunds Universitets Matematiska Seminarium 
{Comm. Sem. Mathem. Universite de Lund) Tome suppl. dedie a M. Riesz, 
Lund (1952) 73-81. 

[6] B.J. West and V. Seshadri, Linear systems with Levy fluctuations, Physica A 
113 (1982) 203-216. 

[7] R. Gorenflo and F. Mainardi, Random walk models for space-fractional diffusion 
processes, Fractional Calculus and Applied Analysis 1 (1998) 167-191. 



11 



[8] R. Gorenflo and F. Mainardi, Approximation of Levy-Feller diffusion processes 
by random walk models, J. for Analysis and its Applications, submitted. 

[9] R. Gorenflo and F. Mainardi, A collection of random walks for Levy-Feller 
diffusion processes, in preparation. 

[10] B.B. Mandelbrot, Fractals and Scaling in Finance (Springer- Verlag, New York, 
1997). 

[11] R.N. Mantegna and H.E. Stanley, Physics investigation of financial markets, 
in: F. Mallamace and H.E. Stanley, eds., Proceedings International School of 
Physics "E. Fermi", Course CXXXIV (10 Press, Amsterdam, 1997) 473-489. 

[12] J.-P. Bouchaud and M. Potters, Theorie des Risques Financiers (Comm. 
Energie Atomique, Alea Saclay, Paris, 1997). 

[13] E. Lukacs, Characteristic Functions (Griffin, London, 1960). 

[14] J.E. Gillis and G.H. Weiss, Expected number of distinct sites visited by a 
random walk with an infinite variance, J. Math. Phys. 11 (1970) 1307-1312. 

[15] B.V. Gncdenko and A.N. Kolmogorov, Limit Distributions for Sums of 
Independent Random Variables (Addison- Wesley, Cambridge, Mass., 1954). 

[16] J.M. Chambers, C.L. Mallows and B.W. Stuck, A method for simulating stable 
random variables, J. Amer. Stat. Assoc. 71 (1976) 340-344. 

[17] R. Bartels, Generating non- normal stable variates using limit theorems 
properties, J. Stat. Comp. Simulation 7 (1978) 199-212. 

[18] R.N. Mantegna, Fast, accurate algorithm for numerical simulation of Levy 
stable stochastic processes, Phys. Rev. E 49 (1994) 4677-4683. 

[19] A. Janicki and A. Weron, Simulation and Chaotic Behavior of a-Stable 
Stochastic Processes (Marcel Dekker, New York, 1994). 

[20] G. Samorodnitsky and M.S. Taqqu, Stable non-Gaussian Random Processes 
(Chapman & Hall, New York, 1994). 

[21] E.W. Montroll and B.J. West, On an enriched collection of stochastic processes, 
in: E.W. Montroll and J. Leibowitz, eds., Fluctuation Phenomena (North- 
Holland, Amsterdam, 1979), pp. 61-175. 

[22] E.W. Montroll and M.F. Shlesinger, On the wonderful world of random walks, 
in: J. Leibowitz and E.W. Montroll, eds., Nonequilibrium Phenomena II: from 
Stochastics to Hydrodynamics (North-Holland, Amsterdam, 1984), pp. 1-121. 

[23] J. Klafter, M. F. Shlesinger and G. Zumofen, Beyond Brownian motion, Physics 
Today 49 No 2 (1996) 33-39. 

[24] C. Tsallis, Levy Distributions, Physics World 10 No 7 (1997) 42-45. 

[25] F. Mainardi, P. Paradisi and R. Gorenflo, Probability distributions generated by 
fractional diffusion equations, in: J. Kertesz and I. Kondor, eds., Econophysics: 
an Emerging Science (Kluwer, Dordrecht, 1998) in press, pp. 39. 



12 



[26] A. Carpinteri and F. Mainardi, eds., Fractals and Fractional Calculus in 
Continuum Mechanics (Springer Verlag, Wien, 1997). 

[27] G.M. Zaslavsky, Physics of Chaos in Hamiltonian Systems (Imperial College 
Press and World Scientific, London and Singapore, 1998). 

[28] R. Hilfer, ed., Applications of Fractional Calculus in Physics (World Scientific, 
Singapore, 1998). 



13 



